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1.0  Introduction 


Considerable  progress  has  been  made  in  the  understanding  of  magnetospheric 
processes  in  the  last  25  years.  During  this  time  much  of  the  effort  has  been 
focused  on  understanding  processes  operating  in  the  tail  of  the  magnetosphere 
and  near  the  magnetospheric  boundaries.  The  inner  magnetosphere  has  not 
been  extensively  investigated  in  the  last  10  to  20  years.  The  Combined  Release 
and  Radiation  Effects  Satellite  (CRRES)  Program  is  designed  to  make  a 
substantial  contribution  to  understanding  this  region  of  space. 

In  this  effort  McDonnell  Douglas  Space  Systems  Company  (MDSSC)  will 
introduce  novel  new  modeling  approaches  and  will  satisfy  one  of  main  goals  of 
the  CRRES  mission,  the  development  of  new  static  and  dynamic  radiation  belt 
models.  Such  models  are  an  important  tool  for  engineers  for  the  design  of 
systems  that  can  survive  in  the  space  environment.  The  present  radiation  belt 
models  that  are  now  used  by  both  the  scientific  and  engineering  communities 
are  the  Vette  radiation  models  developed  by  the  National  Space  Sciences  Data 
Center  (NSSDC)  in  the  late  1960’s  and  early  70’s.  The  data  sets  which  were 
used  to  develop  the  Vette  models  were  acquired  by  instmments  that  are  quite 
primitive  compared  to  today's  state-of-the-art  instruments.  Nevertheless  these 
older  models  have  served  the  community  well. 

The  present  NSSDC  developed  models  are  organized  in  B,  L  space,  a 
coordinate  system  developed  by  C.  E.  Mcliwain  in  1961  (Mclilwain,  1961).  This 
coordinate  system  has  been  virtually  unchanged  since  that  time.  Improvements 
have  come  only  in  the  form  of  improved  computational  techniques.  Although  the 
B,  L  system  has  proven  useful  in  the  inner  zone,  its  use  in  the  outer  zone  has  not 
been  as  successful. 
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This  CRRES  analysis  effort  will  include  the  development  of  new  tools  for  the 
organization  of  the  data.  This  effort  will  develop  novel  new  techniques  for 
organizing  charged  particle  data  in  the  inner  and  outer  zone.  In  the  inner  zone 
we  will  develop  a  model  that  not  only  takes  into  account  the  effect  of  the 
magnetic  field  in  organizing  the  charged  particles  but  also  the  effect  of  the  solar 
cycle  dependent  atmosphere  in  shaping  the  low  altitude  region  of  the  inner 
radiation  belt.  In  the  outer  zone  we  will  provide  a  coordinate  system  that  can 
correctly  represent  adiabatic  changes  in  the  radiation  belt  and  fully  take  into 
account  drift  shell  splitting  and  yet  represent  the  entire  outer  zone  in  terms  of 
only  two  parameters,  an  equivalent  L  and  an  effective  equatorial  pitch  angle. 
Once  the  new  tools  developed  under  this  effort  are  implemented  and  a  best  fit  is 
made  to  the  CRRES  data,  a  high  quality  radiation  belt  model  will  be  created  that 
will  be  valid  for  ail  epoch  within  a  solar  cycle  and  for  all  magnetic  conditions  of 
the  magnetosphere.  In  the  inner  zone,  the  new  model  will  permit  calculation  of 
the  fluxes  during  any  part  of  the  solar  cycle.  In  the  inner  and  outer  zone  the 
model  will  help  separate  adiabatic  changes  from  non  adiabatic  variations,  allow 
the  calculation  of  particle  fluxes  for  all  states  of  magnetospheric  compression, 
and  will  help  theorists  explain  many  of  the  observed  changes  within  the 
magnetosphere. 

One  of  the  most  important  contributions  of  the  CRRES  satellite  is  its 
demonstration  that  the  inner  magnetosphere  is  considerably  more  dynamic  than 
previously  suspected.  Many  satellites  have  investigated  the  region  in  the  tail  and 
the  region  near  synchronous  orbit.  The  large  variations  in  this  region  are  much 
appreciated  and  have  been  the  focus  of  numerous  investigations.  On  the  other 
hand  the  inner  magnetosphere  (below  L  =  3  or  4)  has  been  ignored  because  of 
its  apparent  stability  and  lack  of  interesting  'events'.  The  March  1991  event,  in 
which  CRRES  showed  a  very  large  permanent  change  to  the  inner  radiation  belt 
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came  as  a  complete  surprise  to  many.  It  should  not  have,  except  for  the  lack  of 
a  detailed  study  in  this  region  of  space. 

The  newly  discovered  dynamics  of  the  inner  radiation  belt  place  an  additional 
burden  on  the  modeling  techniques.  In  order  to  develop  realistic  models  of  the 
radiation  environment,  it  is  imperative  that  the  processes  that  create  the 
enhancements  are  understood.  Since  events  such  as  the  March  1991  event  are 
very  rare  and  isolated,  a  single  event  can  be  studied  in  considerable  detail 
without  contamination  by  other  events.  Furthermore,  the  size,  of  the  particle 
energization,  the  size  of  the  magnetic  and  electric  fields  make  the  identification 
of  the  various  signatures  more  precise.  This  event  is  the  first  well-documented 
event  in  this  class.  It  is  now  apparent  that  it  is  not  the  first  and  most  certainly  not 
the  last  event  of  its  class. 

Magnetic  and  electric  field  models  as  well  as  particle  motion  models  are 
important  tools  in  studying  the  dynamics  of  the  inner  magnetosphere.  Some  of 
the  tools  necessary  to  study  this  event  were  developed  many  years  ago  but  were 
never  needed  or  used.  In  particular,  a  vector  potential  model  of  the 
magnetospheric  magnetic  field  developed  almost  15  years  ago  was  more  of  a 
curiosity  since  the  magnetic  field  derived  from  it  was  not  as  accurate  as  a  more 
direct  model  of  the  magnetic  field.  This  vector  potential  model,  does  however, 
become  an  important  tool  in  understanding  the  electric  field  during  the  March 
1991  event. 

The  probability  of  understanding  this  event  is  much  greater  than  understanding 
the  sub-storm  problem.  Although,  magnetic  field  variations  are  very  large, 
several  hundred  nanotesla,  and  the  electric  field  exceeds  several  hundred 
milliVolts/meter,  the  total  magnetic  field  is  over  a  thousand  nanotesla  in  this 
region  of  space  and  thus  the  coordinate  system  as  defined  by  the  magnetic  field 
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is  relatively  stable.  This  is  very  important  simplification  in  our  task  of 
understanding  the  magnetosphere. 
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2.0  An  Electromagnetic  Model  of  the  Magnetosphere  During  the 
March  Event 

The  next  sections  describe  the  beginnings  of  a  modeling  effort  that  we  believe 
has  a  very  good  probability  of  understanding  a  very  important  magnetospheric 
process  using  basic  principles  of  physics. 

2.1  A  Vector  Potential  Model  of  the  Magnetopause  Current  System 

The  first  step  in  the  development  of  the  1977  quiet-time  Olson  Pfitzer  magnetic 
model  was  the  definition  of  the  current  systems  (Olson  and  Pfitzer,  1 977).  The 
magnetopause  current  system  was  represented  by  a  set  of  triangular  current 
sheets.  The  ring  and  tail  currents  were  represented  by  a  set  of  wire  loops.  An 
integration  was  then  performed  over  the  current  systems  to  calculate  both  the 
magnetic  field  and  the  magnetic  vector  potential  from  the  three  current  systems 
at  a  large  number  of  points  within  the  magnetosphere.  The  integration  effort 
produced  a  set  of  files  containing  the  magnetic  field  and  the  vector  potential  at  a 
large  number  of  locations  for  various  dipole  tilt  values.  Orthonormal  function 
fitting  techniques  were  then  used  to  derive  functional  forms  for  both  the  vector 
potential  and  the  magnetic  field.  The  curl  of  the  vector  potential  model  was  then 
compared  to  calculated  magnetic  field  values.  The  curl  of  the  vector  potential 
gave  a  magnetic  field  that  was  a  reasonable  representation  of  the  calculated 
magnetic  field  values.  However,  the  functions  determined  directly  from  the 
magnetic  field  values  produced  a  magnetic  field  model  that  was  considerably 
superior  to  the  magnetic  field  model  using  the  curl  of  the  vector  potential  model. 
We  thus  decided  to  set  aside  the  vector  potential  model  and  published  the 
functional  form  of  the  direct  fit  magnetic  field  model.  This  became  the  much 
used  1977  tilt  dependent  model  for  quiet  time.  This  model  is  valid  to  13  Re  and 
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has  been  used  extensively  by  many  researchers.  It  has  been  extensively 
validated  and  found  to  have  considerable  precision  during  quiet  times. 

A  somewhat  later  effort  used  the  above  described  vector  potential  model  to 
study  the  magnetic  and  electric  field  variations  from  the  wobbling  dipole.  The 
induction  electric  field  calculated  from  the  tilt  dependent  vector  potential  was 
used  to  study  possible  particle  accelerations  due  to  the  wobbling  dipole.  This 
effect  proved  to  be  quite  small,  although  non-trivial.  The  vector  potential  model 
was  thus  relegated  to  dark  deep  desk  drawer. 

It  is  this  vector  model,  developed  almost  1 5  years  ago,  that  was  the  starting  point 
for  the  present  analysis.  One  should  remember  that  the  current  systems  that  are 
used  for  its  definition  are  the  same  current  systems  that  are  used  for  the  highly 
successful  1977  magnetic  field  model.  The  original  vector  potential  model 
developed  in  1977  contained  the  effects  of  all  of  the  current  systems  and  thus 
the  functions  used  in  the  fit  are  unnecessarily  complicated  for  this  work.  For  the 
study  of  the  acceleration  due  to  the  changing  magnetopause  currents,  one 
should  use  only  the  magnetopause  currents.  It  was  possible  to  easily  separate 
the  coefficients  for  the  magnetopause  currents.  Redefining  the  functional  form 
to  remove  terms  that  help  fit  the  ring  current  would  have  been  very  labor 
intensive.  Although  simpler  functions  would  increase  computational  efficiency 
and  accuracy,  the  precision  required  for  this  initial  study  did  not  justify  this 
additional  work.  Thus  the  vector  potential  model  used  for  this  study  consist  of  a 
set  of  polynomials  and  polynomials  times  an  exponential  that  has  virtually  the 
same  form  as  the  1 977  magnetic  field  model.  The  coefficients  for  the  model  are 
the  coefficients  derived  from  the  magnetopause  currents.  As  discussed  above, 
the  accuracy  of  this  vector  potential  model  was  validated  in  1977  when  the  curl 
of  the  vector  potential  was  compared  point  for  point  against  the  magnetic  field 
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values  calculated  from  the  current  systems  and  the  total  calculated  magnetic 
field  was  compared  to  the  delta  B  contours  of  Sugiura  (Sugiura,  et.  al.,  1971). 
The  vector  potential  model  listed  in  Appendix  A  is  thus  a  high  fidelity  model  of 
the  vector  potential  developed  from  the  1 977  tilt  dependent  model  and 
specifically  modified  for  this  effort  to  include  only  the  magnetopause  currents. 

2.2  Extending  the  Vector  Potential  Model  to  Disturbed  Times 

The  vector  potential  model  developed  from  the  1977  tilt  dependent  current 
system  is  only  valid  during  quiet  time.  This  model  was  extended  for  use  during 
disturbed  times  using  techniques  developed  during  the  Consolidated  Data 
Analysis  Workshops  (CDAW).  Extensive  work  with  the  various  (CDAW)  data 
sets  validated  a  method  of  extending  a  quiet  time  model  to  disturbed  times.  This 
method  has  been  shown  to  work  particularly  well  for  scaling  the  magnetopause 
currents  in  response  to  changes  in  the  stand-off  distance.  The  justifications  for 
the  scaling  techniques  are  discussed  in  Olson  and  Pfitzer  1982.  This  same 
scaling  technique  can  be  applied  to  the  vector  potential  model.  Appendix  A 
describes  a  FORTRAN  code  for  a  quiet  time  magnetopause  vector  potential 
model,  a  disturbed  vector  potential,  a  quiet  time  magnetic  field  as  derived  from 
the  curl  of  the  vector  potential  and  a  disturbed  magnetic  field  model  as 
calculated  using  the  curl.  Appendix  A  also  describes  the  algorithms  for 
extending  the  quiet  time  model  to  disturbed  times. 

2.3  The  Induction  Electric  Field 

I  he  great  advantage  of  a  vector  potential  model  is  that  one  can  calculate  the 
Induction  electric  field  directly.  The  induction  electric  field,  is  given  by 
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Note  that  all  work  in  this  analysis  is  performed  in  MKS  units. 

To  calculate  the  time  dependent  electric  field  one  must  have  a  time  dependent 
vector  potential.  The  vector  potential  routines  described  in  Appendix  A  give  the 
vector  potential  as  a  function  of  the  standoff  distance.  Thus  a  time  dependent 
standoff  distance  will  produce  a  time  dependent  vector  potential.  This  vector 
potential  can  then  be  used  to  calculate  a  time  dependent  in'^uction  electric  field. 

2.4  Pressure  Balance  During  the  March  1991  Event 

To  develop  a  model  for  the  magnetic  field  and  for  the  induction  electric  field 
during  the  large  sudden  corr.mencement  of  March  1991 ,  we  need  a  t  me  history 
of  the  solar  wind  pressure  during  the  event.  Unfortunately  very  little  solar  wind 
information  is  available  for  ihis  event.  We  thus  assume  that  before  the  sudden 
commencement,  the  solar  wind  pressure  was  nominal  and  that  the  standoff 
distance  was  at  10.5  Re.  We  also  assume  that  there  was  a  discontinuous 
change  in  the  solar  wind  pressure,  a  step  function  in  the  solar  wind  pressure. 
Travel  times  of  the  solar  wind  plasma  from  sun  suggest  that  the  velocity  of  the 
shock  during  the  March  1991  event  was  on  the  order  of  1450  km/sec.  We 
further  assumed  that  the  magnetopause  was  compressed  to  a  minimum 
distance  of  about  5  Re.  The  magnetooause  cannot  instantly  respond  to  a 
pressuie  discontinuity  and  thus  some  time  required  for  the  standoff  distance  to 
change  from  10.5  to  5  Re.  The  time  dependent  standoff  distance  lunuuon  was 
first  developed  such  that  the  calculated  magnetic  signature  matched  the 
signature  observed  by  the  CRRES  magnetometer.  The  CRRES  magnetometer 
sees  a  positive  dB/dt  lasting  for  approximately  30  seconds  (see  figure  1 ).  Thus, 
the  initial  work  assumed  that  the  standoff  distance  moved  from  10.5  Re  to  5  Re 
in  30  seconds.  The  time  rate  of  change  of  the  standoff  distance  change  was 
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assumed  to  be  nonlinear.  Since  the  magnetic  field  becomes  stronger  as  the  field 
is  compressed,  the  rate  of  change  of  distance  should  slow  down  as  the  boundary 
moves  in.  Several  functional  forms  were  constructed.  The  final  function  was 
selected  on  the  basis  of  correctly  modeling  the  dB/dt  observed  by  the  CRRES 
magnetometer.  The  function  that  was  found  to  give  an  acceptable  result  has  the 
form, 


R,  =  10.5 


(2) 


where  Rj  is  the  standoff  distance  and  t  is  the  time  since  the  start  of  the  event. 
Figure  2  gives  the  time  dependent  standoff  distance  predicted  by  the  above 
equation.  Using  this  standoff  distance  function  for  the  vector  potential  model 
gives  a  dB/dt  at  the  location  of  CRRES  as  shown  in  figure  3.  This  compares 
very  favorably  with  the  CRRES  data  shown  in  figure  1 . 

The  above  determined  of  rate-of-change  of  standoff  distance  has  no  theoretical 
foundation.  George  Siscoe  of  UCLA  (private  communication)  noted  that  during  a 
sudden  change  in  solar  wind  pressure,  the  magnetosphere  must  remain  in 
equilibrium  at  all  times.  That  is,  when  the  solar  wind  pressure  changes,  the 
magnetospheric  boundary  must  move  so  as  to  maintain  pressure  balance  during 
the  motion.  The  velocity  term  in  the  pressure  balance  equation  is  not  the 
velocity  of  the  solar  wind,  but  the  difference  in  the  velocity  of  the  solar  wind  and 
the  velocity  of  the  moving  boundary.  Thus 

R.=98.[pV^]-'"'  (3) 


where  V  is  the  velocity  difference  between  the  solar  wind  speed  and  the  velocity 
of  the  boundary.  The  constant  98,  is  the  constant  that  was  developed  for  our 
dynamic  magnetic  field  models.  This  can  be  rewritten  to  give 
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where  is  the  solar  wind  velocity,  dr /dt  is  the  velocity  of  the  boundary,  y  is 
6371  and  changes  the  boundary  velocity  from  units  or  Rg/sec  to  km/sec,  r  is  the 
instantaneous  position  of  the  boundary.  This  gives  rise  to  the  following  simple 
differential  equation 


dt  = 


_  yr^  dr 


This  can  be  written  as  an  easy  to  solve  integral 
f  r  dr 
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where  a  =  98^  /^p  and  b  =  -  v 
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Evaluating  the  integral  gives 


r=R, 

(8) 

r=10  5 

where  k3=a/b 

This  gives  the  time,  t,  since  the  start  of  the  arrival  of  the  solar  wind  pressure 
change  as  a  function  of  the  instantaneous  location  of  the  standoff  distance  R^. 
Since  no  value  was  available  for  p  the  solar  wind  number  density,  we  used  a 
value  of  30  which  is  consistent  with  a  minimum  standoff  distance  of  5  at  a  solar 


t  =  -^ 
a 


__k_ 

r-10.5  3b 


2  k^-kr  +  r"  VIk 


-10- 


wind  velocity  of  1450  km/sec.  Substituting  values  into  the  above  equation  gives 
the  time  dependent  standoff  distance  as  a  function  of  time  determined  using 
dynamic  pressure  balance.  The  result  of  this  dynamic  pressure  balance  analysis 
is  given  in  figure  4.  We  note  with  interest  that  this  figure  is  very  similar  to  figure 
2.  Figure  2  is  determined  by  attempting  to  fit  the  CRRES  magnetometer 
observations  and  figure  4  attempts  to  use  a  more  theoretical  approach.  Both 
methods  could  be  substantially  improved  if  actual  measurements  become 
available  of  the  real  solar  wind  velocity  and  particle  density  during  this  event. 

This  analysis  is,  however,  consistent  enough  to  allow  us  some  confidence  that 
either  figure  2  or  4  can  be  used  as  the  starting  point  for  the  development  of  a 
time  dependent  vector  potential.  For  this  report  we  have  used  the  much  simpler 
form  (equation  2).  The  time  dependent  vector  potential  is  driven  by  the  time 
dependent  standoff  distance.  This  allows  us  to  investigate  the  induction  electric 
field  created  by  the  change  in  location  and  the  change  in  strength  of  the 
Chapman  Ferraro  currents.  During  this  compression,  the  Chapman  Ferraro 
currents  move  inward  from  10.5  to  5.5  Rq  and  increase  in  strength  by  an  order 
of  magnitude  in  a  time  period  of  approximately  30  seconds. 

2.5  Magnetic  signature  During  the  March  Event 

Figure  3  gives  the  derivative  of  the  magnetic  signature  determined  by  the  time 
dependent  magnetopause  model.  The  integral  of  this  magnetic  field  change,  the 
amplitude  of  actual  delta  B  spike  predicted  by  the  model  at  the  location  of 
CRRES,  is  75  nanotesla.  At  2.7  R©  on  the  noon  meridian,  the  delta  B  spike  is 
predicted  to  be  240  nanotesla.  At  the  surface  of  the  earth,  the  magnetic  field 
spike  should  vary  from  120  nanotesla  at  midnight  to  170  nanotesla  at  local  noon. 
Mid-latitude  magnetometers  may  see  a  magnetic  field  spike  that  will  exceed  this 
magnitude  due  to  currents  induced  in  the  earth.  The  induction  currents  due  to 
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the  conducting  earth  may  increase  the  magnetic  field  signature  as  much  as  60 
percent.  The  exact  enhancements  have  not  yet  been  calculated.  We 
understand  the  size  of  the  increases  to  the  Sq  signatures  that  have  a  period  of 
24  hours.  It  is  not  prudent,  however,  to  apply  the  same  increase  to  a  change 
with  the  period  of  30  seconds.  In  order  to  determine  the  induction  currents  in  the 
surface  of  the  earth  one  must  solve  the  problem  of  a  magnetized  conducting 
sphere  with  finite  conductivity  during  a  30  second  magnetic  field  pulse.  This  is  a 
non-trivial  problem  that  may  require  extensive  analysis. 

2.6  The  Induction  Electric  Field  During  the  March  Event 

The  time  dependent  standoff  distance  given  in  figure  2  was  used  to  calculate  the 
induction  electric  field  during  this  event  at  various  locations  within  the 
magnetosphere.  Figure  5  gives  the  induction  electric  field  at  CRRES  during  the 
30  second  compression  period.  One  notes  the  rapid  rise  in  the  electric  field  to  a 
level  of  approximately  50  mV/m.  The  present  model  only  represents  the  period 
of  active  inward  motion.  At  the  end  of  the  compression  period  the  electric  field 
will  rapidly  decay  to  zero  and  then  actually  reverse  since  the  magnetosphere 
relaxed  somewhat  after  the  initial  compression  (see  figure  1 ).  Figure  6  gives  the 
induction  electric  field  on  the  local  noon  meridian  at  a  distance  of  2.7  Re  from  the 
center  of  the  Earth.  One  notes  that  this  field  increases  to  almost  400  mV/m. 

The  rate  of  change  of  the  induction  electric  field  in  figure  6  differs  from  the  rate  of 
change  seen  in  figure  5.  This  is  due  to  the  fact  that  figure  6  is  on  the  noon 
meridian  and  much  closer  to  the  approaching  currents.  The  electric  field  at  this 
location  is  due  not  only  to  the  increasing  strength  of  the  Chapman-Ferraro 
currents  but  also  to  the  rapidly  approaching  currents.  CRRES  which  is  toward 
the  dark  side  of  the  earth  is  farther  from  the  currents  and  thus  is  most  sensitive 
to  changes  in  the  strength  of  the  Chapman-Ferraro  currents.  Figure  7  is  a 
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snapshot  of  the  electric  field  in  the  equatorial  plane  at  time  t  =  20  seconds.  At 
this  time  the  magnetopause  boundary  is  passing  through  6.0  R©.  The  induction 
electric  field  is  given  every  1 .0  Rg  on  a  rectangular  grid.  The  length  of  the  line  is 
proportional  to  the  strength  of  the  field.  A  line  .5  Re  long  corresponds  to  a  field 
value  of  500  mV/m.  The  direction  of  the  line  gives  the  direction  of  the  induction 
electric  field  vector.  The  dot  at  the  end  of  the  line  points  in  the  positive  field 
direction.  We  note  that  the  field  is  in  most  places  tangent  to  the  azimuthal 
direction  and  that  the  direction  of  the  field  on  the  sunlit  side  is  such  that  both 
electrons  and  protons  will  experience  a  large  gain  in  energy.  There  will  be 
deceleration  near  local  midnight,  but  this  is  small  since  the  field  is  much  smaller 
in  this  region.  Once  the  compression  of  the  magnetic  field  stops  and  the 
magnetosphere  relaxes  the  electric  field  pattern  will  reverse  its  direction.  From 
the  magnetometer  data  in  figure  1 ,  one  can  see  that  only  a  partial  relaxation 
occurs  and  thus  the  deceleration  fields  will  be  weaker.  Furthermore  any  particle 
that  was  near  noon  during  the  start  of  the  acceleration  will  most  likely  be  near 
local  midnight  during  the  deceleration  phase  and  will  thus  be  shielded  from  the 
deceleration  phase.  The  induction  electric  field  is  a  non-conservative  field.  Even 
if  the  acceleration  and  deceleration  fields  were  equal  and  opposite,  some  of  the 
particles  would  experience  substantial  permanent  acceleration.  There  would  of 
course  be  a  class  of  particles  that  would  experience  permanent  deceleration. 

2. 7  The  Parallel  Electric  Field 

Because  of  symmetry,  the  induction  electric  field  from  the  magnetopause 
currents  is  perpendicular  to  the  magnetic  lines  of  force  in  the  magnetic  equatorial 
plane.  At  all  other  locations  there  is  a  component  of  the  induction  electric  field 
that  is  parallel  to  the  field  lines.  Since  the  conductivity  along  field  lines  is  very 
high,  there  will  be  a  very  rapid  redistribution  of  charges  along  the  line  of  force 
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such  that  the  total  electric  field  parallel  to  the  lines  of  force  is  zero.  The  total 
electric  field  is  given  by 

ri\  — 

=  (9) 

ot 

-VO  is  the  scalar  potential  electric  field  is  due  to  charge  separation.  Charges 
will  realign  themselves  such  that  the  parallel  component  of  Ej  is  everywhere 
zero.  This  charge  rearrangement  which  cancels  the  parallel  portion  of  the  total 
field  can  substantially  modify  the  electric  field  perpendicular  to  the  lines  of  force. 

It  is  possible  to  calculate  the  charge  separation  electric  field.  During  our  work 
with  the  induction  electric  field  due  to  wobbling  dipole  we  developed  a  routine 
that  performed  a  line  integral  along  a  line  of  force  into  the  ionosphere.  The  total 
electric  field  along  this  line  of  force  was  required  to  remain  zero  everywhere 
along  the  line  of  force.  This  required  introducing  a  potential  variation  along  the 
line  of  force  so  that  the  gradient  of  the  potential  along  the  line  of  force  would 
everywhere  cancel  the  parallel  component  of  the  induction  electric  field. 

Adjacent  line  integrals  can  then  give  the  gradients  of  the  potential  electric  field 
perpendicular  to  the  lines  of  force.  Since  potentials  are  arbitrary  to  within  a 
constant,  the  analysis  depends  on  the  correct  use  of  the  boundary  condition. 
Since  the  foot  of  the  field  line  is  anchored  in  the  conducting  ionosphere,  the 
ionosphere  becomes  the  physical  boundary  condition.  During  our  wobbling 
dipole  analysis  we  used  either  an  equipotential  ionosphere  or  an  ionospheric 
boundary  condition  that  assumed  that  the  earth  was  a  rotating  conducting 
magnetized  sphere. 

A  similar  analysis  will  be  instructive  in  this  case.  This  part  of  the  electric  field 
analysis  has  not  yet  been  completed.  It  is  the  intention  to  place  a  high  priority  on 
this  analysis.  During  our  wobbling  dipole  analysis,  the  induction  electric  field  was 
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very  small  and  the  difference  between  the  two  boundary  conditions  was  very 
large.  Since  the  driver  induction  electric  field  for  this  event  is  two  orders  of 
magnitude  greater  than  that  of  the  wobbling  dipole  field,  we  expect  less 
sensitivity  to  the  form  of  the  ionospheric  boundary  condition.  We  do,  however, 
expect  a  substantial  change  in  the  overall  electric  field  pattern.  In  many  cases, 
canceling  the  parallel  electric  field  may  substantially  increase  the  perpendicular 
electric  field. 
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3.0  Particle  Acceleration 


Since  the  electric  field  is  approximately  azimuthal  in  the  equatorial  plane  one  can 
make  a  quick  estimate  of  the  amount  of  energy  gain  that  one  can  expect  for 
protons  and  electrons.  The  energy  gain  is  simply  the  E  dl.  Estimating  the  path 
length  from  local  dawn  to  dusk  or  dusk  to  dawn  and  multiplying  by  approximately 
400  mV/m  gives  an  approximate  energy  gain  of  20  MeV  at  R  =2.7  and  30  MeV 
at  R  =  4.  These  are  very  large  numbers  and  suggest  that  the  induction  electric 
field  due  to  the  Chapman  Ferraro  currents  is  very  important  in  understanding 
the  particle  energization  that  CRRES  observed  during  the  March  1991  event. 

3. 1  Lorentz  Force 

The  force  on  a  charged  particle  is  given  by 

F  =  q(E  +  vxB)  (10) 

A  modified  cosmic  ray  trajectory  code  was  used  to  integrate  the  trajectory  of 
protons  using  equation  10.  The  initial  cosmic  ray  code  used  by  many 
Investigators  was  modified  to  step  In  time  instead  of  position.  It  was  modified  to 
include  the  effects  of  the  electric  field  and  energy  conservation  was  removed 
from  the  code.  The  code  can  perform  on  the  order  of  50,000  integration  steps 
before  round  off  errors  begin  to  affect  the  accuracy  of  the  code.  Thus  proton 
motion  during  the  event  can  easily  be  studied.  However,  the  motion  of  electrons 
cannot  easily  be  studied  by  a  trajectory  integration  program.  To  study  electron 
motion  a  guiding  center  code  must  be  used.  For  this  analysis  only  proton 
trajectories  were  studied. 
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3.2  Particle  Motion  During  the  March  Event 


A  Lorentz  force  trajectory  code  is  very  attractive  because  of  its  simplicity. 

Cosmic  ray  codes  have  been  extensively  verified  and  shown  to  be  accurate. 

The  Lorentz  force  equation  includes  all  effects.  For  our  analysis  we  used  the 
electric  field  as  given  by  the  induction  electric  field  due  to  the  changing 
Chapman-Ferraro  currents.  The  electric  field  is  calculated  from  the  time 
dependent  vector  potential.  Similarly  the  magnetic  field  consists  of  a  dipole  field 
plus  the  time  dependent  magnetic  field  calculated  from  the  curl  of  the  same  time 
dependent  vector  potential.  The  magnetic  and  electric  field  codes  are  described 
in  Appendix  A.  An  overview  of  the  Lorentz  force  integration  code  is  given  in 
Appendix  B.  The  summary  of  the  Lorentz  force  code  is  presented  in  Appendix  B 
for  reasons  of  completeness.  The  listing  of  the  code  will  allow  the  user  to  easily 
verify  the  results  of  the  analysis  presented  in  this  document. 

To  perform  the  actual  analysis,  the  trajectory  code  was  reversed  and  a  negative 
proton  was  integrated  backward  in  time.  This  allows  us  to  begin  the  study  at  the 
point  of  observation  and  determine  the  initial  starting  point  of  the  particle  at  time 
zero.  One  of  the  important  milestones  in  the  analysis  of  the  March  event  will  be 
the  development  a  model  of  particle  motion  and  acceleration  that  can  show  the 
origin  of  the  particles  seen  by  the  CRRES  spacecraft.  Are  the  particles 
accelerated  from  the  local  population  or  are  they  accelerated  inward  from  the 
cosmic  ray  flux  in  the  outer  zone? 

Figure  8  shows  a  sample  trajectory  calculation.  The  particle  was  started  at 
approximately  3  R©  and  3  hours  local  time  with  an  energy  of  50  MeV  and  a  pitch 
angle  of  90  degrees.  The  particle  was  started  at  time  t  =  30  seconds.  At  time 
t  =  0  seconds  the  proton  was  at  a  local  time  of  21  hours  and  had  an  energy  of  10 
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MeV  and  was  at  a  radial  distance  of  4.5  Re-  Thus  during  the  30  seconds  of 
positive  dB/dt,  this  test  particie  drifted  through  270  degrees.  From  Figure  7  one 
can  see  that  almost  the  entire  drift  period  was  in  a  region  where  the  electric  field 
was  in  a  direction  necessary  for  acceleration.  Many  other  trajectories  were 
analyzed.  The  amount  of  energy  gain  is  strongly  associated  with  the  drift 
velocity.  The  particle  whose  drift  velocity  is  such  that  it  drifts  at  least  1 80 
degrees  in  30  seconds  shows  the  largest  amount  of  energy  gain.  Very  low 
energy  particles  have  very  low  drift  velocities  and  thus  their  E  dI  is  very  small 
since  the  drift  path  length  in  30  seconds  is  short.  These  particles  will  also  be 
decelerated  by  the  relaxation  of  the  boundary  after  the  initial  increase.  The 
particles  with  fast  drift  velocity  will  have  drifted  to  the  night  side  of  the  earth 
where  the  deceleration  effect  is  small,  but  the  slow  particles  will  still  be  on  the 
dayside  and  will  experience  the  full  deceleration. 

Table  1  is  a  summary  of  results.  The  table  lists  starting  and  ending  Ls  as  well  as 
starting  and  ending  energies.  From  this  table  one  can  see  that  energy  gain  is 
larger  at  higher  energies  and  that  the  energy  gain  is  larger  at  larger  distances. 
Furthermore,  the  change  in  L  is  less  for  the  smaller  L  shells.  A  quick 
investigation  of  this  table  shows  that  protons  with  a  final  L  shell  of  2.4  originated 
from  Ls  in  the  vicinity  of  3.0  Re.  These  protons  most  likely  originate  from  the 
existing  trapped  proton  flux.  Protons  with  a  final  L  of  3.0  originate  from  Ls 
greater  than  4.5  and  could  thus  have  their  origin  in  the  solar  proton  flux  present 
in  the  outer  zone  at  this  time. 

The  results  in  Table  1  are  all  for  particles  with  a  pitch  angle  of  90  degrees. 
Several  trajectories  were  run  for  non  90  degree  particles.  It  is  apparent  from 
these  runs  that  the  energy  gain  is  in  the  component  of  energy  perpendicular  to 
the  magnetic  field.  Thus  this  acceleration  mechanism  will  produce  particles  with 
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4.0  Discussion  and  Summary 

This  report  details  the  first  steps  in  developing  a  storm  time  model  of  the 
magnetic  and  electric  fields  during  a  very  large  sudden  compression  of  the 
magnetosphere.  We  believe  that  the  model  of  the  Chapman-Ferraro  currents  for 
this  event  is  a  good  approximation  for  this  event.  Two  separate  techniques,  one 
based  on  reproducing  the  CRRES  magnetic  signature,  and  one  based  on 
maintaining  dynamic  equilibrium  during  compression  give  approximately  the 
same  result.  The  biggest  uncertainty  is  in  the  actual  solar  wind  pressure 
parameters. 

One  fundamental  simplification  has  been  made  in  the  analysis  of  the  Chapman- 
Ferraro  currents.  In  the  analysis  we  assume  that  the  entire  magnetosphere 
responds  as  a  unit.  The  entire  magnetosphere  contracts  at  the  same  time  and 
at  the  same  rate,  and  the  entire  current  system  increases  its  value  uniformly 
over  the  entire  magnetopause.  In  reality,  the  front  of  the  magnetopause  will 
change  first  with  the  more  distant  tail  currents  responding  much  later.  However, 
the  most  important  currents,  the  currents  that  dominate  the  magnetic  signature 
in  the  region  of  interest  (the  inner  magnetosphere),  are  within  a  few  Rg  of  the 
sub  solar  point.  Thus  this  approximation  does  not  introduce  a  very  large  error. 
The  biggest  error  in  the  analysis  to  date  is  that  we  have  assumed  a  vacuum 
magnetosphere.  We  have  not  yet  included  the  effect  of  the  plasma  within  the 
magnetosphere.  An  additional  complication  of  a  non-vacuum  magnetosphere  is 
that  the  electromagnetic  disturbance  created  by  the  changing  Chapman-Ferraro 
currents  will  propagate  with  a  finite  speed,  most  likely  the  Alfven  speed.  This 
may  introduce  delays  in  the  arrival  of  the  disturbance  vector  proportional  to  the 
distance  from  the  sub  solar  point.  This  propagation  delay  should  be  on  the  order 
of  10s  of  seconds  and  thus  should  be  measurable  by  using  timing  marks  at 
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various  magnetometer  locations  on  satellites  and  on  the  ground.  The  rise  time 
of  the  event  is  very  sharp  and  thus  delay  time  measurements  with  an  accuracy 
of  seconds  can  reasonably  be  expected.  Furthermore,  since  CRRES  measures 
both  the  magnetic  and  electric  field  disturbances,  it  should  be  possible  to  easily 
calculate  the  direction  of  the  propagation  vector  of  the  disturbance  pulse. 

This  event  offers  a  unique  opportunity  for  investigating  the  injection  of  particles 
into  the  trapping  region  using  only  Maxwell's  equations.  The  injection  is  in  a 
region  where  the  magnetic  field  is  strong  and  the  magnetic  coordinate  system 
not  extensively  distorted.  CRRES  is  in  the  very  center  of  the  injection  region. 
The  magnetopause  vector  potential  model  provides  an  excellent  first  look  insight 
into  the  energy  that  is  available  from  the  changing  magnetopause  current 
system.  The  first  look  energization  shown  in  Table  1  show  that  more  than 
enough  energy  is  available  and  that  the  injection  mechanism  is  the  induction 
electric  field.  One  must  now  solve  the  very  complex  details  of  the  problem.  Two 
of  the  most  important  details  are  the  fact  that  the  magnetosphere  does  contain  a 
plasma  and  that  the  electromagnetic  wave  does  have  a  finite  propagation 
velocity.  One  of  the  first  priorities  during  the  next  year  of  this  effort  will  be  the 
reevaluation  of  the  electric  field  using  the  E  -  B  =  0  condition.  We  expect  that 
including  this  effect  in  the  electric  field  model  will  change  the  electric  field  pattern 
and  that  the  electric  field  perpendicular  to  the  magnetic  lines  of  force  will  be 
increased. 

The  CRRES  data  set  in  conjunction  with  various  other  data  sets  are  of  sufficient 
quality  to  greatly  assist  in  the  development  of  a  complete  model  of  this  unique 
event.  The  final  model  of  the  March  1991  event  will  be  an  important  milestone  in 
the  understanding  of  magnetospheric  dynamics. 
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5.0  Appendix  A 

Tnis  appendix  describes  and  lists  a  series  of  subroutines  that  describe  the  magnetic 
and  electric  fields  during  the  March  1991  event. 

5. 1  Subroutine  AXYZ 

This  subroutine  is  the  basic  sector  potential  subroutine.  The  routine  was  developed  in 
1977  along  with  the  1977  Olson  Pfitzertilt  dependent  magnetic  field  model.  The 
routine  calculates  the  magnetic  vector  potential  in  units  of  nanotesla-Re  everywhere 
inside  the  magnetosphere  and  inside  of  a  sphere  of  radius  13  Re.  The  routine  is  tilt 
dependent.  The  routine  is  a  series  of  polynomials  plus  polynomials  times  an 
exponential.  The  complexity  of  the  function  is  such  that  if  the  coefficients  of  the  ring 
and  tail  are  added  to  coefficients  describes  in  this  listing,  the  functions  have  sufficient 
fidelity  to  describe  the  detailed  structure  due  to  the  ring  and  tail  current  systems. 

5.1.1  Calling  sequence 

XX(3)  a  3  dimension  input  array  that  specifies  the  position  in  Cartesian  solar  magnetic 
coordinates.  XX(3)  along  the  north  dipole  axis,  XX(1 )  is  perpendicular  to  XX(3) 
and  in  the  plane  containing  XX(3)  and  the  sun-earth  line  and  pointing  in  the 
direction  of  the  sun,  XX(2)  completes  the  right  handed  coordinate  system.  The 
distance  are  given  in  unit  of  Re. 

AT(3)  a  3  dimensioned  array  that  returns  the  vector  components  of  the  magnetic 
vector  potential.  The  units  are  in  nanotesla-Re. 

COMMON/TILTIT/TILT  TILT  is  an  input  variable  that  specifies  the  tilt  of  the  earth's 
dipole  axis.  Zero  tilt  indicates  that  the  dipole  is  perpendicular  to  the  sun-earth 
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line.  Positive  tilt  is  when  the  northern  dipole  is  tipped  toward  the  sun.  This 
value  must  be  set  up  before  a  call  is  made  to  routine  AXYZ. 


5.1.2  Program  Ustina  -  Subroutine  AXYZ 


c 

c 

C 

C 


SUBROUTINE  AXYZ (XX, AT) 

This  routine  caluclates  the  Vector  potential  of  the  magnetopause 
magnetospheric  magnetic  field  during  quiet  conditions  for  any  tilt. 
XX  (3)  is  a  real*8  position  in  earth  radii  in  solar  magnetic  coords 
AT (3)  is  the  real*8  vector  potential  in  nanotesla-Re 
REAL*8  TILT 

REAL* 8  D(44) ,E(64) ,F(44) ,DM(88) ,EM(128) ,FM(88) , 

*AT(3) ,X(10) ,y(10) ,Z(10) ,XX(3) ,TT(4) ,AA(3) , 

*TILTL, XN, YN, ZN, R2, R 
INTEGER*2  ITD (44) , ITE  (64) , ITF  (44) , 

*1, II, K 

COMMON/TILTIT/TILT 


DATA  ITD  /1,2, 1, 
*2, 1,1, 1,2, 1,2,1, 
DATA  ITE  /I, 2, 1, 
*1,1, 1,2, 1,2, 1,2, 
*1, 1,2/ 

DATA  ITF  /2,1,2, 
*1,2,2,2,1,2,1,2, 
DATA  (DM(I),I=1, 
*-.729348268D+00, 
*-.651629108D-01, 
*-.313618130D-02, 

*  .100224170D-04, 
*-.3026819320-03, 

*  .236910507D-05, 

*  .  395551400D-05, 

*  .328361431D-05, 

*  .344430885D-06, 

*  .245904885D-05, 
*-.2128370260-05, 
*44*0./ 

DATA  (EM (I), 1=1, 

*  .  539534465D+01, 

*  .848876852D+00, 
*-.5983766030-01, 

*  .422206531D-04, 

*  .201674318D-02, 
*-.7417480710-05, 

*  .459399443D-01, 

*  .562030059D-03, 
*-.2656416100-05, 

*  .558513796D-06, 

*  .  149384635D-05, 


2, 1,1, 2, 1,1, 1,2, 1,2, 1,1, 2, 2, 1,2, 1,1, 1,1, 2, 1,2, 1,1, 
1,2, 2, 1,2, 1,1,1/ 

2, 1,1, 2, 1,1, 1,2, 1,1, 2, 1,1, 2,  2, 2, 2, 2,  1,2,  1,1, 1,1, 2, 
1,1, 2, 1,1, 1,2, 1,1, 2, 1,1, 2, 2, 2, 2, 2, 1,2, 1,1, 1,1, 2,1, 


1,2, 2, 1,2, 2, 2, 1,2, 
2, 1,1, 2, 1,2, 2, 2/ 
88)  / 

-.126711994D-03,-. 
-.166931408D-04,-. 
-.569972419D-05,  . 
.4474169040-08,  . 
.526556403D-06,  . 
.404562068D-09,-. 
.1964637640-07,-. 
-.3811503200-08,  . 
.1339821190-08,  . 
.9862300490-09,  . 
-.8721991840-08,  . 


1,2,2, 1, 1,2, 1,2, 2, 2, 2, 1,2, 1,2,2, 


2502562450-02, 

4290008330-03, 

1015812730-02, 

1571650050-04, 

1399323880-03, 

1749082880-02, 

1720101050-04, 

1887425050-05, 

2457456570-03, 

1834092400-04, 

1721697130-04, 


-.9297340000-06, 

-.2741239690-07, 

.2108758570-05, 

.3714356080-07, 

.1639432760-06, 

.1110436130-05, 

.1043338570-06, 

.7960652640-08, 

.2203851170-06, 

.2651296010-07, 

.2940762300-08, 


128)  / 

-.9172467290-03,- 
.7103295020-04, 
-.1651772820-06,- 
. 1731705870-07, 
.8699260960-05,- 
-.9197579090-08,- 
.1630723910-04, 

- .2630838650-06, - 
-.1733641070-09,- 
-.3555052590-09,  - 
-  .  1841426100-08, - 


.2022122270-02, 

.3050087230-02, 

.3349583320-02, 

.1609667760-03, 

.1704243380-02,  ■ 

.426006943D-05,■ 

.4459576850-03, 

.4002415680-04,  ■ 

.1115868490-03, 

.1734413250-04, 

.5655507020-03,  ■ 


. 1296506340-05, 
. 1205514600-05, 
.4722135530-05, 
.2080889960-08, 
.2801639070-05, 
.5338883390-07, 
.5119261300-08, 
.2280380440-06, 
. 4873506900-07, 
. 1391979380-08, 
. 1271316360-05, 
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*-.3778264570-05, -.24 1832155D-07,  . 3950292000-04, -. 1507960740-06, 

*  .1418634270-04,  . 1575477700-07, -. 8947417770-05,  .4414003720-07, 

*  .869472 6570-05, -.1701708360-07, -.12315 16500-06, -.67 5272 34 50-10, 

*-.5495662500-04, -. 1344172660-07,  .5383396840-05,  .1193980670-07, 

*-. 1058664290-03, -.2244207710-06, -.2989115210-05, -. 4323032360-09, 
*64*0  .  / 

OATA  (FM(I) , 1=1, 88) / 

*  .4520638590-02,  .3093212350-06,  .6581182670-01,  .9225436910-05, 

*  .2219302720-03, -.2693767990-06,  . 4343661160-02, -. 1694772960-05, 

*  .2642680490-03, -.4 965390620-07, -.3830523880-04, -.8187459940-08, 
*-.13 9054 9300-03, -.9027530360-07, -.12204 42540-05,  .3059951140-09, 

*  .2994056010-05, -.8239802520-08, -.1370819260-06,  .6867451070-09, 

*-.4494079540-05,  . 5162724890-08, -. 1088071560-03, -. 2251097120-07, 

*-.3873139660-03, -.1814613770-06, -.257 81 4 4 680-05,  .3621727350-09, 

*  .107887 52 90-05, -.81780101 90-10,  . 3272707560-04 ,-. 2137684900-06, 

*  .1823 9515 90-04, -.4 676910710-07, -.42572 98210-05,  .4575647280-08, 
*-.5074171820-04,  . 2344726820-07, -. 8288387280-06, -. 2224328420-09, 

*  .1184875950-06,  .2595589680-10,  .4450233790-06,  .3457922920-09, 
*44*0./ 

OATA  TILTL/99.O+0/ 

IF(TILT.EQ.TILTL)GO  TO  20 

TILTL=TILT 

TT(1)=1. 

TT(2)=TILT 
TT(3)=TILT*TILT 
TT (4) =TT (3) *TILT 
DO  10  1=1, 64 
II=(I-1) *2+1 
K=ITD(I) 

if (i.le.44)D(I)=(10.*DM(II) ) *TT(K)+(10.*DM(II+1) ) *TT(K+2) 
K=ITE(I) 

E (I)=(10. *EM(II) ) *TT (K) +(10.*EM(II+1) ) *TT (K+2) 

K=ITF (I) 

10  if (i . le. 44) F (I) = (10. *FM(II) ) *TT (K) + (10 . *FM (II+l) ) *TT (K+2) 

20  CONTINUE 
XN=XX ( 1 ) 

YN=XX(2) 

ZN=XX(3) 

R2=XN**2+YN**2+ZN**2 
R=SQRT (R2) 

DO  1  1=1,7 
X(I)=XN 
Y ( I ) =YN 
Z (I) =ZN 
XN=XN*XX (1) 

YN=YN*XX(2) 

ZN=ZN*XX(3) 

1  CONTINUE 
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*+D(21)*Y(  5)+D(22)*X(  4)*Y(  1) 

AA(1)  =AA{1)  +  (0 . 0  +D(23)*Y(  1)+D(24)*Y(  1)*Z(  1) 

*+D(25)*X(  1)*Y(  1)+D(26)*X(  l)*y(  1)*Z(  1)+D(27)*Y(  1)*Z(  2)+D(28) 
**Y(  3)+D(29)*Y(  3)*Z(  1)+D(30)*y(  3)*Z(  2)+D(31)*X(  l)*y(  1)*Z(  2) 
*+D(32)*X(  l)*y(  3)+D(33)*X(  l)*y(  3)*Z(  1)+D{34)*X(  2)*y(  1)+D(35) 

**X(  2)*Y(  1)*Z(  1)+D(36)*X(  2)*y(  1)*Z(  2)+D(37)*X(  2)*y(  3)+D(38) 

**y(  1)*Z(  3)+D(39)*X(  l)*y(  1)*Z(  3)+D(40)*X(  3)*y(  1)+D(41)*X(  3) 

**Y(  1)*Z(  1)+D(42)*y(  1)*Z(  4)+D(43)*y(  5)+D(44)*X(  4)*y{  1))*EXP 

*(  -.06*R2) 

AA(2)=+E(  1)+E(  2)*Z(  1)+E{  3)*X(  1)+E(  4)*X(  1)*Z{  1)+E(  5)*Z(  2) 

*+E(  6)*Y(  2)+E(  7)*y(  2)*Z(  1)+E(  8)*y(  2)*Z(  2)+E(  9)*X(  1)*Z(  2) 

*+E(10)*X(  l)»y(  2)+E(ll)*X(  l)*y(  2)*Z(  1)+E(12)*X(  1)*Y(  2)*Z(  2) 

*+E(13)*X(  2)+E(14)*X(  2)*Z(  1)+E(15)*X(  2)*Z(  2)+E(16)*X(  2)*Y(  2) 

*+E(17)*X(  2)*y(  2)*Z(  1)+E(18)*Z(  3)+E(19)*y(  2)*Z(  3)+E(20)*X(  1) 
**Z(  3)+E(21)*X(  2)*Z(  3)+E(22)*X(  3)+E(23)*X(  3)*Z(  1)+E(24)*X(  3) 
**Z(  2)+E(25)*X(  3)*y(  2)+E{26)*Z(  4)+E(27)*Y(  4)+E(28)*y(  4)*Z(  1) 
*+E(29)*X(  1)*Z(  4)+E(30)*X(  l)*y(  4>+E(31)*X(  4)+E(32)*X(  4)*Z(  1) 
AA(2)=AA(2) 

*+(0.0  +E(33)+E(34)*Z(  1)+E(35)*X(  1)+E(36)*X(  1)*Z(  1)+E(37)*Z(  2) 
*+E(38)*Y{  2)+E(39)*y(  2)*2(  1)+E(40)*Y(  2) *Z  (  2)+E(41)*X(  1)*Z(  2) 

*+E(42)*X(  l)*y(  2)+E(43)*X(  l)*y(  2)*Z(  1)+E(44)*X(  l)*y(  2)*Z(  2) 

*+E(45)*X(  2)+E(46)*X(  2)*Z(  1)+E(47)*X(  2)*Z(  2)+E(48)*X(  2)*y(  2) 

*+E(49)*X(  2)*y(  2)*Z(  1)+E(50)*Z(  3)+E(51)*Y(  2)*Z(  3)+E(52)*X{  1) 
**Z(  3)+E(53)*X(  2)*Z(  3)+E(54)*X(  3)+E(55)*X(  3)*Z(  1)+E(56)*X(  3) 
**Z(  2)+E(57)*X(  3)*Y(  2)+E(58)*Z(  4)+E(59)*y(  4)+E(60)*Y(  4)*2(  1) 
*+E(61)*X(  1)*Z(  4)+E(62)*X(  l)*y(  4)+E(63)*X(  4)+E(64)*X(  4)*Z(  1) 
*)*EXP  (  -.06*R2) 

AA(3)=+F(  l)*y(  1)+F(  2)*y<  1)*Z(  1>+F(  3)*X(  l)*y(  1)+F(  4)*X(  1) 

**y(  1)*Z(  1)+F(  5)*Y<  1)*Z(  2)+F(  6)*y(  3)+F(  7)*Y(  3)*Z(  1)+F(  8) 

**y(  3)*Z(  2)+F(  9)*X(  l)*y(  1)*Z(  2)+F(10)*X(  1)*Y(  3)+F(11)*X(  1) 

**y(  3)*Z(  1)+F(12)*X(  2)*Y(  1)+F(13)*X(  2)*Y(  1)*Z(  1)+F(14)*X(  2) 
**Y(  1)*Z(  2)+F(15)*X(  2)*Y(  3)+F(16)*y(  1)*Z(  3)+F(17)*X{  1)*Y(  1) 
**Z(  3)+F(18)*X(  3)*y(  1)+F(19)*X(  3)*y(  1)*Z(  1)+F(20)*Y(  1)*Z(  4) 
*+F(21)*y(  5)+F{22)*X(  4)*y(  1) 

AA(3)=AA(3)+(0.0  +F(23)*Y(  1)+F(24)*y(  1)*Z(  1) 

*+F(25)*X(  l)*y(  1)+F(26)*X(  l)*y(  1)*Z(  1)+F(27)*y(  1)*Z(  2)+F(28) 
**y(  3)+F(29)*y(  3)*Z(  1)+F(30)*y(  3)*Z(  2)+F(31)*X(  l)*y(  1)*2(  2) 
*+F(32)*X(  l)*y(  3)+F(33)*X(  l)*y(  3)*Z(  1)+F(34)*X(  2)*y(  1)+F(35) 

**X(  2)*y(  1)*Z(  1)+F(36)*X(  2)*y(  1)*Z(  2)+F(37)*X(  2)*y(  3)+F{38) 

**y(  1)*Z(  3)+F(39)*X(  l)*y(  1)*Z(  3)+F(40)*X{  3)*y(  1)+F(41)*X(  3) 

**Y(  1)*Z(  1)+F(42)*Y(  1)*Z(  4)+F(43)*y(  5)+F(44)*X(  4)*y(  1) ) *EXP 
*(  -.06*R2) 

AT ( 1 ) =AA ( 1 ) 

AT(2)=AA(2) 

AT{3)=AA(3) 

RETURN 

END 
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5.2  Subroutine  AXYZDN 


Subroutine  AXYZDN  is  the  disturbed  condition  vector  potential  model.  The  routine 
calls  subroutine  SMAG  which  provides  the  correct  scaling  parameters  for  this  routine. 
The  subroutine  provides  the  correct  disturbed  time  vector  potential  at  time  t  providing 
subroutine  SMAG  provides  the  correct  scaling  parameters.  The  scaling  parameters 
that  are  used  are  STRMAG,  the  strength  of  the  magnetopause  currents,  and  SCL  the 
size  scaling  parameter.  These  parameters  are  discussed  in  section  5.6 

5.2.1  Calling  Sequence 

XX{3)  a  3  dimension  input  array  that  specifies  the  position  in  Cartesian  solar  magnetic 
coordinates.  XX(3)  along  the  north  dipole  axis,  XX(1 )  is  perpendicular  to  XX(3) 
and  in  the  plane  containing  XX(3)  and  the  sun-earth  line  and  pointing  in  the 
direction  of  the  sun,  XX(2)  completes  the  right  handed  coordinate  system.  The 
distance  are  given  in  unit  of  Re. 

AT(3)  a  3  dimensioned  array  that  returns  the  vector  components  of  the  disturbed 
condition  magnetic  vector  potential.  The  units  are  in  nanotesla-Re 

T  The  time  during  the  event.  The  time  must  be  in  units  of  seconds.  The  time,  T, 
is  passed  through  to  SMAG,  where  SMAG  must  use  it  to  determine  the  scaling 
parameters. 

COMMON/TILTIT/TILT  TILT  is  also  an  input  variable.  It  specifies  the  tilt  of  the  earth’s 
dipole  axis.  Zero  tilt  indicates  that  the  dipole  is  perpendicular  to  the  sun-earth 
line.  Positive  tilt  is  when  the  northern  dipole  is  tipped  toward  the  sun.  This 
value  must  be  set  up  before  a  call  is  made  to  routine  AXYZDN. 

5.2.2  Subroutine  Listing  -  AXYZDN 

subroutine  axyzdyn (x, a, t) 
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C  This  is  the  disturbed  condition  Vector  potential.  I  uses  the  scaling 
C  algorithms  developed  in  1982  (JGR  Aug.  82  p5943) 

C  It  requires  that  subroutine  SMAG  return  the  magnetopause  current 
C  and  magnetopause  scale  size  as  a  function  of  time 
C  X(3)  is  the  position  in  Re  in  solar  magnetic  coord 
C  A (3)  is  the  Vector  potential  in  nanotesla-Re 
C  T  is  the  time  in  seconds 

REAL*8  a (3) ,x(3) , xx(3) ,aa (3) , STRMAG, SCL, T 
INTEGER* 2  I 

call  smag(strmag, scl, t) 
do  110  i=l,3 
110  XX (i) =x (i) *scl 

call  axyz(xx,aa) 
do  120  i=l,3 
120  a (i) =aa ( i) *strmag 
return 
end 
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5.3  Subroutine  CURLA 


This  subroutine  calculates  the  curl  of  the  quiet  time  vector  potential.  It  produces  the  a 
quiet  time  magnetic  field  model,  that  except  for  the  precision  in  the  fit  and  numerical 
derivatives  will  be  very  similar  to  the  Olson  Pfitzer  1977  tilt  dependent  model. 


5.3.1  Calling  Sequence 

XX(3)  a  3  dimension  input  array  that  specifies  the  position  in  Cartesian  solar  magnetic 
coordinates.  XX{3)  along  the  north  dipole  axis,  XX{1 )  is  perpendicular  to  XX(3) 
and  in  the  plane  containing  XX(3)  and  the  sun-earth  line  and  pointing  in  the 
direction  of  the  sun,  XX(2)  completes  the  right  handed  coordinate  system.  The 
distance  are  given  in  unit  of  Re. 

BBB(3)  a  3  dimensioned  array  that  returns  the  vector  components  of  the  quiet  time 
magnetic  field.  The  units  are  in  nanotesla. 

COMMON/TILTIT/TILT  TILT  is  also  an  input  variable,  it  specifies  the  tilt  of  the  earth’s 
dipole  axis.  Zero  tilt  indicates  that  the  dipole  is  perpendicular  to  the  sun-earth 
line.  Positive  tilt  is  when  the  northern  dipole  is  tipped  toward  the  sun.  This 
value  must  be  set  up  before  a  call  is  made  to  routine  CURLA. 


5.3.2  Program  Listing  -  Subroutine  CURLA 

SUBROUTINE  CURLA (XX, BBB) 

C  This  subroutine  calculates  the  numerical  CURL  of  the  Vector 
C  potential  and  thus  calculates  the  quiet  time  magnetic  field. 
C  It  calls  AXYZ  and  thus  returns  the  value  of  B  in  nanotesla 
C  XX  is  the  Real*8  value  of  the  position  in  Earth  radii 
C  BBB  is  the  Real*8  value  of  the  magnetic  field  in  nanotesla 
C  DEL  is  a  step  size  parameter  for  the  numercial  CURL 
REAL*8  X(3) ,B{3) ,BB(3,3),BBB(3),xx(3) 

*,DEL 

INTEGER*2  I,J,K 
DATA  DEL/0.0001/ 
do  1  i=l,3 
1  x(i)=xx(i) 

CALL  AXYZ(X,B) 
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DO  10  1=1,3 

X(I)=X(I)+DEL 

CALL  AXYZ  (X.BBd,  I)  ) 

10  X(I)=X(I)-DEL 
DO  20  1=1,3 
J=I+2 

J=J- (J-1) /3*3 
K=I  +  1 

K=K- (K-1) /3*3 

20  BBB ( I ) = (BB ( J, K) -B ( J) -BB (K, J) +B (K) ) /DEL 

RETURN 
END 
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5.4  Subroutine  DYNB 


This  subroutine  calculates  the  disturbed  of  dynamic  magnetic  field  values.  It  produces 

the  a  disturbed  time  magnetic  field  model  using  the  curl  of  the  vector  potential. 

5.4.1  Calling  Sequence 

XX(3)  a  3  dimension  input  array  that  specifies  the  position  in  Cartesian  solar  magnetic 
coordinates.  XX(3)  along  the  north  dipole  axis,  XX(1 )  is  perpendicular  to  XX(3) 
and  in  the  plane  containing  xx(3)  and  the  sun-earth  line  and  pointing  in  the 
direction  of  the  sun,  XX(2)  completes  the  right  handed  coordinate  system.  The 
distance  are  given  in  unit  of  Re. 

B{3)  a  3  dimensioned  array  that  returns  the  vector  components  of  the  disturbed  time 
magnetic  field.  The  units  are  in  nanotesla. 

BMAGreturns  the  magnitude  of  the  disturbed  time  magnetic  field  in  units  of  nanotesla 

T  an  input  variable  that  gives  the  time  during  the  event.  This  time,  T,  must  be  in 
units  of  seconds.  The  time,  T,  is  passed  through  to  SMAG,  where  SMAG  must 
use  the  time  to  determine  the  scaling  parameters. 

COMMON/TILTIT/TILT  TILT  is  also  an  input  variable.  It  specifies  the  tilt  of  the  earth’s 
dipole  axis.  Zero  tilt  indicates  that  the  dipole  is  perpendicular  to  the  sun-earth 
line.  Positive  tilt  is  when  the  northern  dipole  is  tipped  toward  the  sun.  This 
value  must  be  set  up  before  a  call  is  made  to  routine  JYNB. 


5>4.2_Proaram  Listing  -  Subroutine  DYNB 

subroutine  dynb(x,b,bniag,  t) 

C  This  is  the  disturbed  condition  magnetospheric  magnetic  field  model 
C  It  determines  the  magnetopause  magnetic  field  as  a  function  of  time 
C  It  requires  that  subroutine  SMAG  determines  the  magnetospheric 
C  current  strength  and  magnetospheric  scaling  parameter  as  a  function 
C  oft ime 

C  X(3)  is  the  position  in  Re  in  solar  magnetic  coord 
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C  B(3)  is  the  mangetic  field  in  nanotesla 
C  BMAG  is  the  mangetic  field  magnitude 
C  T  is  the  time  in  seconds 

REAL*8  x(3) ,xx(3) ,b(3) , bb (3) , BMAG, T, STRMAG, SCL 
INTEGER*2  I 

call  smag (strmag, scl, t) 
do  210  i=l,3 
210  XX (i) =x (i) *scl 

call  curlA(xx,bb) 
do  220  i=l,3 
220  b (i) =bb (i) *strmag 

bmag=dsqrt (b(l) **2+b(2) **2+b<3) **2) 

return 

end 
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5.5  Subroutine  E FIELD 


This  subroutine  calculates  the  induction  electric  field  from  the  changing  Chapman 

Ferraro  currents.  It  produces  an  induction  electric  field  as  a  function  of  time  when  the 

Chapman-Ferraro  currents  are  changing  in  response  to  a  changing  solar  wind 

pressure. 

5.5.1  Calling  Seou^^nce 

XX(3)  a  3  dimension  input  array  that  specifies  the  position  in  Cartesian  solar  magnetic 
coordinates.  XX(3)  along  the  north  dipole  axis,  XX(1)  is  perpendicular  to  XX(3) 
and  in  the  plane  containing  XX{3)  and  the  sun-earth  line  and  pointing  in  the 
direction  of  the  sun,  XX(2)  completes  tho  right  handed  coordinate  system.  The 
distance  are  given  in  units  of  meters. 

E(3)  a  3  dimensioned  array  that  returns  the  vector  components  of  the  disturbed  time 
magnetic  field.  The  units  are  in  Volts/meter. 

EMAGreturns  the  magnitude  of  the  disturbed  time  magnetic  field  in  units  of  Volts/meter 

T  an  input  variable  that  gives  the  time  during  the  event.  The  time  must  be  in  units 
of  seconds.  The  time,  T,  is  passed  through  to  SMAG,  where  SMAG  must  use  it 
to  determine  the  scaling  parameters. 

COMMON/TILTIT/TILT  TILT  is  a  variable  used  by  the  vector  potential  program.  This 
routine  sets  the  value  of  TILT  to  0.  The  tilt  is  the  tilt  of  the  earth’s  dipole  axis. 
Zero  tilt  indicates  that  the  dipole  is  perpendicular  to  the  sun-earth  line.  Positive 
tilt  is  when  the  northern  dipole  is  tipped  toward  the  sun.  This  value  must  be  set 
up  before  a  call  is  made  to  routine  EFIELD. 

5.5.2  Program  Listing  --  Subroutine  EFIELD 

subroutine  ef ield (xx, E, emag, t ) 

C  This  routine  calculates  the  induction  electric  field.  It  mu.st  e 


C  in  MKS  units. 

C  XX (3)  is  the  position  is  entered  in  meters  (solar  magnetic  coords) 

C  E(3)  returns  the  vector  induction  magnetic  field  in  Volts/meter 
C  E  =  negative  of  the  time  derivative  of  the  vector  potential 
C  Emag  is  the  magnitude  of  the  induction  electric  field. 

C  T  is  the  time  in  seconds.  Subroutine  SMAG  must  be  properly  set  up 
C  to  give  the  magnetospheric  boundary  parameters  as  a  function  of  the 
C  time  in  seconds 

REAL*8  X(3),XX(3),al(3),a2(3),E(3) 

REAL*8  EMAG, T, TILT, T1,T2,DELT, CON 

INTEGER*2  I 

COMMON/TILTIT/TILT 

data  con/6. 371D-3/,delt/0. 0001/ 

do  10  i=l,3 

10  X  (i) =xx ( i) /6 . 37 lD+6 
tilt=0 
tl=t 

t2=tl+delt 

call  axyzdyn (x, al, tl) 
call  axyzdyn (X, a2, t2) 

E(l)=-(a2(l)-al(l) )/delt*con 
E  (2) =- {a2 (2) -al (2) ) /delt*con 
E (3) =- (a2 (3) -al (3) ) /delt*con 
emag=dsqrt (E (1) **2+E (2) **2+E (3) **2) 
return 
END 
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5.6  Subroutine  SMAG 


This  routine  is  the  time  dependent  driver  routine  that  must  be  modified  by  the  user  to 
give  the  time  dependent  magnetopause  scaling  factors.  The  routine  must  calculate  the 
strength  of  the  magnetopause  currents,  and  the  scale  size  of  the  magnetosphere  as  a 
function  of  the  time,  t.  Time  must  be  in  units  of  seconds.  The  scale  factor  SCL  is 
given  by 

SCL  =  10.5/Rs 

The  magnetopause  current  strength  factor  STRMAG  is  given  by 


STRMAG  = 


Where  Rs  is  the  standoff  distance 


5.6.1  Calling  Sequence 

T  is  the  time  in  units  of  seconds.  It  is  passed  through  by  the  various  field  routines. 
This  routine  converts  time  to  the  time  dependent  scale  factors  and 
magnetopause  current  strength  factor. 

STRMAG  returns  the  strength  of  the  magnetopause  current  at  time  t 

SCL  is  the  time  dependent  scale  factor  that  scales  the  positions  with  respect  to  the 
size  of  the  magnetopause. 


S,6,Z  Pmrm  listing  -  SubrQuiin&.SMAG. 


subroutine  smag (strmag, scl, t) 

C  This  subroutine  gives  the  standoff  distance  and  scaling  parameters 
C  as  a  funtion  of  the  time  in  seconds 
REAL*  8  STRMAG, SCL, T, STDOFF 
C 

C  This  is  an  example  of  a  linear  change  in  the  standoff  distance  at 
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o  o 


.3  Re/second 
goto  800 

STDOFF=10.5-.3*T 
SCL=10.5/STDOFF 
STRMAG=SCL**3 
return 
C 

C  This  is  an  example  of  a  different  form  for  the  change  in  magnetopause 
C  configuration.  It  tries  to  match  the  CRRES  dB/dt  observation  for 
C  rev  587  during  the  first  30  seconds  of  the  event. 

C 

800  if  (t.le.O)then 
strmag=l 
else 

strmag=l+9* ( (t/30) **1.3) 
endi  f 

scl=strmag**0 .3333333 
stdof f=10 . 5/scl 
end 
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6.0  Appendix  B  -  Lorentz  Force  Integration  Program 

This  section  briefly  describes  a  Runge  Kutta  trajectory  integration  program  It  is 
included  in  this  report  for  completeness  since  it  was  used  to  verify  proton  acceleration. 
It  will  permit  any  user  to  reproduce  the  values  calculated  in  this  document.  The 
Runge-Kutta  techniques  used  in  this  program  are  derived  from  a  cosmic  ray  cut-off 
code  that  is  over  20  years  old.  The  code  which  was  initially  written  to  integrate  the 
path  as  a  function  on  position.  The  code  used  variable  size  distance  steps  and  was 
written  to  hold  the  energy  of  a  particle  fixed.  The  code  was  modified  to  include  the 
elect; ic  field  and  the  distance  stepping  algorithm  was  changed  to  a  time  step  code. 

The  code  has  variable  step  size.  The  step  size  depends  on  the  Larmor  radius  of  the 
particle  trajectory  and  the  drift  rate  of  the  electric  field  force.  The  code  uses  MKS  units 
throughout. 

Since  the  program  was  not  developed  for  general  use,  no  attempt  is  made  here  to 
describe  each  of  the  routines  in  detail.  The  routines  contain  a  substantial  number  of 
comment  cards.  This  along  with  the  overall  simplicity  of  the  code  should  permit  most 
user  to  successfully  reproduce  the  work  in  this  document. 

6. 1  Program  Listing 


PROGRAM  TRAJCHK 

C  This  is  a  driver  program  that  sets  up  a  call  to  the  trajectory  program 
c 

reai*8  X ( 6) , v, t ilt , xmag, t 

real *4  r, xlong, th,ph,w,anl,an2,an 

integer*2  i 

common/t ilt it /t ilt 

C  Set  tilt  angle  to  zero  and  request  starting  coords.  Enter  distance  in  Re 
c  Angle  from  noon  is  +  toward  +y  or  dusk,  pitch  is  angle  with  repeat  to  z=0 
c  plane,  azimuth  is  +  toward  +x 


print  *, 'Enter  R, Angle  from  noon, Azimuth, Pitch, W 
1005  read  * , R, xlong, th, ph, W 

c  Tape  7  writes  a  file  of  the  results,  this  copy  is  set  up  to  inegrate 
c  backward  in  time.  To  change  to  forward  a  sign  must  be  changed  in  3 
c  locations  these  are  flagged  by  CS$$$$$$$$$$$$$$$ 
write  (7,*)'  Backward' 
write (7,  *) r, xlong, th, ph, w 
anl=th*3. 14159/180. 
an2=ph*3 . 14159/180 
call  veloc(W,v) 
print  *,w,v 
an=xlong*3 .14159/180. 

C  Set  up  intial  posistion  X(l)  thru  X(3)  hold  particle  position  in  meters 
c  X(4)  thru  x(6)  hold  particles  intitial  velocity  in  meters/second 
C 

X(l)=r*6.371e+6*cos(an) 

X(2) =r*6 . 371E+6*sin (an) 

X(3)=0 

X (4) =cos (anl) *sin (an2) 

X(5)=sin(anl) *sin (an2) 

X ( 6) =cos  (an2) 

xmag=dsqrt (x(4) **2+x(5) **2+x(6) **2) 
do  5  i=4,6 

5  X (i) =x (i) /xmag*V 

t=30. 

CALL  TRAJPRO(X,t) 

END 

SUBROUTINE  VELOC(W,V) 

C  Determine  the  initial  velocity  of  the  proton  given  its  energy  in  MeV 
c  W  is  the  energy  in  MeV 
c  V  is  the  velocity  in  meter/second 
REAL*8  V,V2C2 

V2C2=W* (W+2*931) / (W+931) **2 
V=3 . 0E+8*DSQRT (V2C2) 

RETURN 

END 

SUBROUTINE  TRAJPRO (X, tt ) 

C  Calculate  one  complete  particle  trajectory  stating  at  position  xx 
C  and  time  t.  XX  is  in  meters  and  time  is  in  seconds 
integer*4  n, number 
real*8  X ( 6) , S (6) , RR ( 6) , Q (6) , t, tt 

rea 1*8  dxdt , xx, bb, bmag, e, emag, eb, RTPF, P2 9, OP 7 , DS, DT, DV, CON, 
*P5DS,P29DS,OP7DS,DIST 
integer*2  numb, i 

real*4  RRR, angle, en, xxx, xxy, xxz, xsv 

COMMON/SAD/DXDT (6) ,XX(6) ,BB(3) ,BMAG,E (3) ,EMAG,EB, t 

common/plotit/xsv (3, 5000) , number 

t=tt 

EB=1  . 

numb=300 

nutnber=0 

N=0 
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DIST=0 
CON=.02 
DO  10  1=1,6 
XX(I)=X(I) 

10  Q(I)=0 

RTPF=DSQRT(.5D0) 

P29=l .DO-RTPF 
OP7=l .DO+RTPF 
CALL  EBFORCE 
EB=EMAG/BMAG 
C 

C  Main  inegration  loop  —  first  set  up  variable  step  size 
C  This  uses  Gill's  method  of  Runge-Kutta 
50  DS=CON*6.57E-8/BMAG 

DV=0.03*SCRT(XX(4) **2+XX(5) **2) 

if  (eb.ne.O)  then 

DT=1 . 05E-8/EMAG*DV 

DS=AMIN1 (DS,DT) 

endif 

P5DS  =  .5*DS 
P29DS=P29*DS 
OP7DS=OP7*DS 
C 

C  GILL'S  NUMERICAL  INTEGRATION  ROUTINE 
C 

DO  60  1=1,6 

SID  =  P5DS*DXDT<I) 

RR(n=S(I)  -Q(I) 

XX(I)=XX{I) +RR(I) 

60  Q(I)=Q(I)+3.*RR(I)-S(I) 

CALL  EBFORCE 

DO  61  1=1, 6 

S(I)  =  P29DS*DXDT(I) 

RR(I)=S(I) -P29*Q(I) 

XX(I)=XX(I) +RR(I) 

61  Q(I)=Q(I)+3.*RR(I)-S<I) 

CALL  EBFORCE 

DO  62  1=1,6 

S(I)  =  OP7DS*DXDT(I) 

RR(T)=S(I) -OP7*Q(I) 

XX(I)=XX(I)+RR(I) 

62  Q(I)=Q(I)+3.*RR(I)-S(I) 

CALL  EBFORCE 

DO  63  1=1,6 

S(I)  =  P5DS*DXDT(I) 

RR(I)=(S(I)-Q(I))/3. 

XX(I)=XX(I)+RR(I) 

63  Q(I)=Q(I)+3.*RR(I)-S(I) 

N=N+1 

DIST  =  DIST  +  DS 

C$$$$$S$$$$  to  change  to  forward  in  time  change  sign  to  +ds 
t=t-ds 

CALL  EBFORCE 

RRR=DSQRT(XX(1) **2+XX(2) **2+xx<3) **2) /6.371E+6 
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c 

C  Every  so  often  print  information  on  progress  to  the  screen 
IF(N/100*100.EQ.N)  then 
angle=datan2 (xx(2) ,xx(l) ) *180-/3.14159 
EN= (XX (4) *XX(4) +XX(5) *XX(5) +xx(6) **2) /I . 912E+14 
WRITE (*, 90) N, t , rrr, emag, angle, en 
endif 

90  FORMAT(I7,f7.2,2fl0.5,2fl0.3) 

C 

C  Write  stuff  to  a  file  every  NUMB  steps 
IF (N/NUMB*NUMB . EQ . N) THEN 
XXX=XX(1) /6.371E+6 
XXY=XX(2) /6.371E+6 
xxz=xx (3) /6 . 371e+6 

EN=(XX(4) *XX(4)+XX(5) *XX ( 5 ) +xx ( 6) * *2 ) /I . 912E+1 4 
angle=datan2 (xx(2) ,xx(l) ) *180-/3.14159 
WRITE  (7,101)  t, RRR, XXX, XXY,xxz,EN, emag, angle 
ENDIF 

101  FORMAT (8f 11 . 3) 

100  FORMAT(I5,4E15.8,/,20X,3E15.8) 

C  check  to  see  if  we  are  still  in  magnetosphere  and  that  time  is 
C  still  valid. 

C  This  exit  condition  must  be  changed  for  forward  in  time  integration 
C$$$$$$$$$$S$$S$$ 

IF  (XXX. GT. -13.0  -and.  rrr. It. 5 
* .and.t .gt .0)  GOTO  50 
RETURN 
END 

SUBROUTINE  EBFORCE 

C  This  routine  calculates  the  Lorentz  Force  on  a  particle 
C  DXDT (1) , (2) , (3)  are  the  drivative  of  the  position 
C  DXDT(4), (5), (6)  are  the  derivative  of  the  velocity 
C  all  equations  use  MKS  units 
C 

REAL* 8  DXDT, XX, BB, BMAG, E, EMAG, EB, t, dcon 
COMMON/SAD/DXDT(6) ,XX(6) , BB (3) ,BMAG, E (3) ,EMAG,EB,t 
0$$$$$$$$$$  to  change  to  forward  in  time  make  constant  + 

DATA  DCON/-9.5D+7/ 

DXDT (1) =XX(4) 

DXDT(2)=XX(5) 

DXDT (3) =XX (6) 

CALL  BFIELD(XX,BB,BMAG,t) 
call  efield(xx,e,emag, t) 

C$$$$S$$$$S$$  to  change  to  forward  in  time  make  change  to  +E  in  next  3 
20  DXDT (4) =DCON* (-E (1) +XX(5) *BB(3)-XX(6)  *BB(2)  ) 

DXDT (5) =DCON* (-E (2) +XX(6) *BB(1) -XX(4) *BB(3) ) 

DXDT (6) =DCON* (-E (3) +XX(4) *BB(2) -XX(5) *BB(1) ) 

RETURN 

END 

SUBROUTINE  BFIELD (X, B, BMAG, t ) 

C  This  routine  combines  a  dipole  field  with  a  dynamic  external  field 
C  It  returns  the  magnetic  field  in  MKS  units  (tesla) 


lines 
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C  XX  the  input  position  is  in  MKS  unites  (meters) 
REAL*8  X,  B,  BMAG,  A,  R,  R2,  t ,  CC,  xx,  btemp,  bb 
integer*2  i 

DIMENSION  X(l) ,B(1) ,xx(3),bb(3) 

DATA  A/-8.1D+15/,cc/1.0D-9/ 
do  300  i=l,3 

300  XX (i) =x ( i) /6 . 371e+6 

call  dynb (XX, bb, btemp, t) 

R2=X(1) **2+X(2) **2+X(3) **2 
R=DSQRT(R2) 

C  Convert  external  field  to  tesla  and  add  dipole 
b(l)=3.d+0*x(l) *x(3) *a/ (r2*r2*r) +bb(l) *cc 
b(2)=3 .d+0*x(2) *x(3) *a/ (r2*r2*r) +bb(2) *cc 
B(3)=(3.d+0*x(3) **2-r2) *A/ (r2»R2*R) +bb (3) *cc 
BMAG=DSQRT(B(1) **2+B (2 ) **2+B (3) **2) 

RETURN 

END 


-40- 


6.0  Figures  and  Tables 


CRRES  magnotometer  data.  Partial  derivative  with  respect  to  time.  Note:  field  ri.ses  rapidly  with  time 
»r  30  seconds,  then  rate  of  change  reverses  direction  but  not  completely. 
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Figure  3.  Time  rate  of  change  of  the  magnetospheric  magnetic  field  at  the  location  of  CRRES  This  field  change  is  calcu¬ 
lated  using  the  standoff  distance  change  shown  in  fiure  2 
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maintained  during  boundary  position  change 
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Figure  5.  Induction  electric  field  at  the  location  of  CRRES  Since  CRRES  is  far  from  the  boundary,  the  primary  contibu 
tion  to  this  change  is  the  change  in  strength  of  the  Chapman-Ferraro  currents 


Figure  6.  Induction  electric  field  on  the  sun-earth  line  at  a  distance  of  2  7  Re  from  the  earth  Change  is  due  to  increase 
in  Chapman-Ferraro  currents  and  to  the  reduced  distance  from  the  observation  point  to  the  currents 
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Figure  7  Electic  field  in  the  magnetic  equatorial  plane  Time  is  t  =  20  seconds 
Standoff  distance  is  at  6  Re.  Electric  field  is  evaluated  on  a  1  0  Re  grid  A  vector  1  0 
Re  long  corresponds  to  a  field  of  0  5  V/m.  Dot  at  the  head  of  the  vector  points  in  the 
direction  of  the  field,  +X  is  down  and  toward  the  sun 
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Figure  8.  Sample  proton  trajectory  This  50  MeV  proton  is  started  at  3  hr  local  time  at 
3  0  Re  The  proton  trajectory  is  integrated  backward  for  30  seconds  The  proton  origi¬ 
nates  at  a  local  time  of  about  21  hrs  local  time  with  an  energy  of  10  MeV  The  proton 
gam  an  energy  of  40  MeV  +X  is  down  and  toward  the  sun 
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Table  1 
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